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ABSTRACT 

We present a new method, the Coupled Escape Probabihty (CEP), for exact calculation of line 
emission from multi-level systems, solving only algebraic equations for the level populations. 
The CEP formulation of the classical two-level problem is a set of linear equations, and we 
uncover an exact analytic expression for the emission from two-level optically thick sources 
that holds as long as they are in the "effectively thin" regime. In comparative study of a number 
of standard problems, the CEP method outperformed the leading line transfer methods by 
substantial margins. 

The algebraic equations employed by our new method are already incorporated in nu- 
merous codes based on the escape probability approximation. All that is required for an exact 
solution with these existing codes is to augment the expression for the escape probability with 
simple zone-coupling terms. As an application, we find that standard escape probability cal- 
culations generally produce the correct cooling emission by the Cll 158 /im line but not by 
the ^P Hnes of O I. 

Key words: radiative transfer — line: formation — methods; numerical — ISM; lines and 
bands 



1 INTRODUCTION 

Much of the information about astronomical sources comes from 
spectral lines, requiring reliable analysis of multi-level line emis- 
sion. Most current methods for exact solutions involve acceler- 
ated A-iteration (ALI) techniques^ in which the radiation inten- 
sity is obtained from the repeated action of an operator designated 
A on the source function (e.g., Rybicki 1991, Hubeny 1992). The 
ALI method utilizing short characteristics with parabolic interpola- 
tion of the source function (hereafter SCP; Olson, Auer & Buchler 
1986) is a standard against which the efficiency of other line trans- 
fer techniques can be measured. 

Because of the complexity and computational demands of ex- 
act methods, many simulation codes that attempt to implement as 
many realistic physical ingredients as possible are altogether by- 
passing solution of the radiative transfer equation, employing in- 
stead the escape probability technique. In this approach only the 
level populations are considered, calculated from rate equations 
that include photon escape factors which are meant to account ap- 
proximately for the effects of radiative transfer (see Dumont et al 
2003 for a recent discussion and comparison with ALI calcula- 
tions). Since this approach is founded on a plausibility assumption 
right from the start, its results amount to an uncontrolled approx- 
imation without any means for internal error estimates. Neverthe- 



less, this inherent shortcoming is often tolerated because of the sim- 
plicity and usefulness of the escape probability approach. 

We present here a new exact method, the Coupled Escape 
Probability (CEP), that retains all the advantages of the naive es- 
cape probability approach. In this new technique the source is di- 
vided into zones, and formal level population equations that are 
fully consistent with radiative transfer are derived rigorously from 
first principles. Different zones are coupled through terms resem- 
bling standard escape probability expressions, resulting in a set of 
level population equations with non-linear coefficients. Solution of 
this set of coupled algebraic equations produces level populations 
that are self-consistent with the line radiation they generate. Any 
desired level of accuracy can be achieved by increasing the number 
of zones. 

We introduce our new method in §2. In §3 we study the 2-level 
model in both a semi-infinite atmosphere and finite slabs, present- 
ing results and comparison with SCP calculations. The new CEP 
method attains the exact solutions, outperforming the SCP method 
by substantial margins. We present the equations for multi-level 
systems in §4 and include as an example an application to the '^P 
system of 01. Section 5 contains a discussion that, among other 
things, covers various technical details. 



^ As noted by Trujillo Bueno & Fabiani Bendicho (1995), the ALI method 
is based on the Jacobi iteration (Jacobi 1845). 



2 THE NEW TECHNIQUE 

Consider the transfer of a line with frequency uq. The dimension- 
less line profile is ^{x), normalized so that J ^{x)dx = 1, where 
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X = — vq)I t^VD is the dimensionless frequency shift from line 
center and IS.vo is the Doppler width. We consider here only the 
case of <I>(a;) that does not change its shape throughout the source 
and frequency independent line source function 5. These assump- 
tions are adopted to simplify the presentation. They do not reflect 
inherent limitations of our new method. 



■b = cos" \l 




X = X, 



x = 



2.1 Radiative Transfer 

For the geometry we adopt a plane-parallel slab whose physical 
properties vary only perpendicular to the surface. Optical depth at 
frequency v along a path orthogonal to the surface is Ti, — r<I>(a;), 
and r can be used as a coordinate that uniquely specifies locations 
in the slab (see figure 0. The optical depth along a ray slanted at 
6 = cos~^/i from normal is r„{fi) — T${x)/fj,, and the intensity 
along the ray obeys the radiative transfer equation 



(1) 



The equation for the flux = 2-k J Ii^fidfi is obtained from in- 
tegration over angles. The overall line flux F — J F^dv obeys at 
every position in the slab 



dF{r] 

dr 
where 



= A-kAud[S{t) - J(r)] 



dQ. 



/ ^ / h{r,^,)^{x)dx 



(2) 



(3) 



is the intensity averaged over both angles and line profile. Denoting 
by Tt the overall optical thickness and accounting for the emission 
from both faces of the slab, the line contribution to the cooling rate 
per unit area is 



A = F(rt) - F(0) EE A-kAvdJ 



(4) 



The line cooling factor j is introduced for convenience when Avd 
is constant in the slab. Integrating equation|2|over r yields 



j= / S{r)p{T)dr 
Jo 

Here we introduced 

S(t) 



(5) 



(6) 



a quantity that has been called the net radiative bracket (Athay & 
Skumanich 1971). From the formal solution of the radiative transfer 
equation, 



p(r) = 1 



when there is no external radiation entering the slab. 



2.2 Level Populations 

Denote by nfc(T), with fc = 1, 2, the populations per sub-state of a 
given transition at position r; that is, nt ~ Nk/gk where gt is the 
level degeneracy and is the overall level population. Then the 
line source function is 



S = 



A2 



B21 ni — n2 



(8) 



zone # 
z 



optical depth 

tz = tt 

— ti 



"^2 
To = 



Figure 1. Top: Sketch of the slab geometry for the radiative transfer prob- 
lem. Bottom: The partition of the slab into zones (see J23). 



where A and B are the Einstein coefficients of the transition. The 
populations are obtained from steady-state rate equations of the 
form ^ Rij = 0. The term corresponding to exchanges between 
the transition levels, separated by E21 — hvo, is 

R21 = -^2in2 - B21 J(n2 - ni) ~ C2i(n2 - me-^^i/'^'^) (9) 

where C is the collision rate; exchanges with other levels have sim- 
ilar form and are listed in 341 



2.3 Solution 

The common approach of exact solution methods is to handle ra- 
diative transfer and the level population distribution as two dis- 
tinct problems, coupled through the results each of them gives. 
The problem is initialized with populations (and the correspond- 
ing source functions) obtained in some limiting case, e.g., ther- 
mal equilibrium. With these populations, radiative transfer (eq.0 
is solved for the intensity to determine J (eq. |3j, which is then 
plugged into the rate terms (eq.|9j to determine new populations, 
and so on. However, from eqs.|6|and|8| the rate term can be written 
as 



R21 — —A2in2p — C21 {n2 — n\e 



-E2i/kT 



(10) 



showing that the only radiative quantity actually needed for the cal- 
culation of level populations at every position is the net radiative 
bracket p(r); given this factor we could compute the level popu- 
lations that are consistent with the radiation they produce without 
solving for the intensity. And as is evident from equations |7| and 
[S] the factor p(r) itself can be computed from the level popula- 
tions, again without solving for the intensity. Therefore, inserting 
p(t) from equation into the rate terms (eg. I10> produces level 
population equations that properly account for all the effects of ra- 
diative transfer without actually calculating the intensity itself, the 
radiative transfer equation has been incorporated through its formal 
solution in equation^ 

A numerical solution of the resulting level population equa- 
tions requires a spatial grid, partitioning the source into zones such 
that all properties can be considered uniform within each zone. The 
degree of actual deviations from uniformity, and the accuracy of the 
solution, can be controlled by decreasing each zone size through 
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finer divisions with an increasing number of zones. FigureQshows 
the slab partitioning into z zones. The z-th zone, i = 1 . . . z, oc- 
cupies the range n-i < t ^ n, with to = and Tz = rt. The 
optical depth between any pair of zone boundaries is 

T''' = \n-rj\ (11) 

so that the optical thickness of the i-th zone is r'''~^. The tem- 
perature and collision rates are constant in the zone, and the corre- 
sponding rate term for its (constant) level populations is 



i?2 



-A2in2P -G2i(n2-nie 'j 



(12) 



where the superscript i is used as a zone label. The factor p(r) 
varies in the zone and has been replaced by a constant p' that 
should adequately represent its value there, for example = 
I b(T"i) or = p (i[ri + Ti_i])). There are no set 

rules for this replacement other than it must obey p{'''i) when 



0. We choose for the zone average 



p{T)dT 



(13) 



and this choice proved to be very successful in our numerical cal- 
culations. From eq.|7| calculation of p* requires an integration over 
the entire slab, which can be broken into a sum of integrals over 
the zones. In each term of the sum, the zone source function can be 
pulled out of the r-integration so that 



P 



dt 



X 



(14) 



The remaining integrals can be expressed in terms of common func- 
tions. Consider for example 



2T-i,i-l 

dr f dt 



$^dx I e-l^-*l*/^^, (15) 



the contribution of zone i itself to p\ It is straightforward to show 
that/3' = /3(r''^-^), where 



I3{t) = - dt ${x)dx dfie 



(16) 



This function was first introduced by Capriotti (1965); it is the 
probability for photon escape from a slab of thickness r, averaged 
over the photon direction, frequency and position in the slab. The 
contribution of zone j 7^ i to the remaining sum can be handled 
similarly, and the final expression for the coefficient p' is 



where 



T-i.i-l^S'. 

j = l 



1/ '-^ 



+ a 



(17) 



(18) 



and where q''-' — t^'-' i3{t^'-'). The quantity q''-' obeys a*'-' = a-' 



and a''' = 0, therefore M'^ = M^' and M'' 



\^ The first 



Since /3* = M" /t^'^ ^ , the first term could be incorporated into the 
sum in eg. llVl as the j = i term. 



term in the expression for p' is the average probability for photon 
escape from zone i, reproducing one of the common variants of 
the escape probability method in which the whole slab is treated as 
a single zone (e.g., Krolik & McKee 1978). The subsequent sum 
describes the effect on the level populations in zone i of radiation 
produced in all other zones. Each term in the sum has a simple in- 
terpretation in terms of the probability that photons generated else- 
where in the slab traverse every other zone and get absorbed in 
zone i, where their effect on the level populations is similar to that 
of radiation external to the slab (see appendix IaI. 

Inserting the coefficients p' from eq. 1171 into the rate terms 
feQ. ll2t in every zone produces a set of non-linear algebraic equa- 
tions for the unknown level populations nj.. The procedure was 
outlined here only for the diffuse radiation of a single transition; 
we describe the extension to multi-levels in !|4|and the inclusion of 
external radiation in appendix]^ Solution of these equations yields 
the full solution of the line transfer problem by considering only 
level populations'^"^ the computed populations are self-consistent 
with their internally generated radiation even though the radiative 
transfer equation is not handled at all. Once the populations are 
found, radiative quantities can be calculated in a straightforward 
manner from summations over the zones. The emerging intensity 
at direction /i is 

z 

I^irt,!^) = ^ (e-^"'*/^ - e-"'-'*/A') S\ (19) 
1=1 

The flux density emerging from each face of the slab obeys 

z 

z 

-F40) = 27r^ [£3(r'"''°-l>) -S3(r''''$)] (20) 

i = l 

where £3 is the third exponential integral (e.g., Abramowitz & Ste- 
gun 1972). The line cooling coefficient is 



= 2/ [ci —a — Q +Q jo 



(21) 



The solution method just described is exact — the discretized 
equations are mathematically identical to the original ones when 
r''*~^ ^ for every i. As is usually the case, the only approx- 
imation in actual numerical calculations is the finite size of the 
discretization, i.e., the finite number of zones. A desired accuracy 
is achieved when, upon further division, the relative change in all 
level populations is smaller than the prescribed tolerance. 



2.4 Numerical Implementation 

The level populations of all zones are described by a set of non- 
linear algebraic equations. The equations are readily solved by the 
Newton method, which utilizes the Jacobian of the set. Since the 
dependence on the unknown variables is explicit in all the rate 
terms, the Jacobian can be computed from analytic expressions. 
The functions /3 (see eq. I16> . a — t/3 and their derivatives are 
conveniently calculated from the representations 



^ Apmzese et al (1980) proposed somewhat similar equations. They based 
their arguments on probabilistic reasoning and did not offer a formal deriva- 
tion. We thank P. Lockett for bringing this to our attention. 
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Figure 2. Plots of the functions /3 (see eg. 1161 and a = t/3 (eg. 1221 . 



Q(r) = 
a{r) = 



^{x)E2[T<i>{x)]dx 



(22) 



where E2 is the second exponential integral. Figure|2|plots /3 and 
Q. While /3 is a monotonically decreasing function, a is mono- 
tonically increasing and its asymptotic behavior when t —> 00 is 
a ~ Vln r. This divergent behavior does not pose any problems 
because a cancels to first order in the thickness of the zones in the 
linear combinations defining (ea. ll8> . Only the second order 
terms, involving the second derivative a", survive. 

Since our aim is to explore the intrinsic accuracy of our new 
method, the integrals in eg. 1221 were computed repeatedly with an 
80 points Gaussian quadrature to ensure that these integrations do 
not compromise the precision of the outcome. The E„ functions 
were evaluated with a rapidly convergent series from Press et al. 
(1986). The integration range was truncated at a; — ±7, which we 
have verified is sufficient in all cases thanks to the rapid decrease 
with X of the integrands. 

In order to test the new method, the radiative transfer prob- 
lem was also solved using the ALI method for comparison. The 
technique is based on a modified A iteration in which the statistical 
equilibrium equations are linearized via the Rybicky & Hummer 
(1992) preconditioning scheme. The method also takes advantage 
of an operator splitting scheme by introducing an approximate op- 
erator A*, the diagonal of the exact operator A in the formal solu- 
tion J = A[S] of the radiative transfer equation. It has been shown 
that the introduction of this operator leads to an optimal balance 
between the convergence rate and computing time per iteration (Ol- 
son, Auer & Buchler 1986; Carlsson 1991). The ALI calculations 
presented in this paper utilize a formal solver based on the short- 
characteristics scheme with parabolic precision (Olson, Auer & 
Buchler 1986), currently considered the method of choice for com- 
plicated line transfer problems (e.g., Kunasz & Auer 1988; Auer, 
Fabiani Bendicho & Trujillo Bueno 1994; van Noort, Hubeny & 
Lanz 2002; Fabiani Bendicho 2003). With this SCP method, equa- 
tionQwas solved for many frequencies and ray inclinations, and the 
mean intensity computed from angular and frequency integrations 



(eq.|3} by numerical quadratures. To ensure the high precision re- 
quired in this comparative study, the angular integration was done 
with a Gaussian quadrature with 24 points in the variable /j.. The 
frequency integrals were done with trapezoidal integration extend- 
ing to 2; = ±4 with 33 frequency points, which we have verified 
yields the desired precision. 

We proceed now to present solutions and comparisons of the 
newly developed CEP method with the SCP method for a num- 
ber of standard problems. In all the examples we employ uni- 
form physical conditions and the Doppler shape for the line pro- 
file, "1> = ■K~^^'^e~^ ; note that the line center optical depth is then 

To = t/^/tV. 



3 2-LEVEL ATOM 

In the two-level problem, the steady-state rate equation R21 = 
(eq.|9) yields the familiar expression for the source function 



S= (l-e)J + eS(r) 



where B is the Planck function and where 

N 



-E2i/kT\ 



(23) 



(24) 



Here is the density of the collision partners and N^^. the stan- 
dard critical density with a slight modification that incorporates the 
Boltzmann factor correction. Replacing J with p (eq.|6}, the equa- 
tion for the source function becomes 



(1 + iip)S = B, 



where 



V 



N' 
If 



(25) 



this result also follows directly from eq.^|with R21 = 0. When the 
2-level problem is formulated with optical depth as the independent 
variable, it is fully characterized by the two input quantities B{T) 
and e (or, equivalently, rj) specified as functions of r. There is no 
need to specify intrinsic properties of the transition such as, for 
example, E21 or A2i. Instead of solving for the population of each 
of the two levels, this single equation for the unknown S provides 
the complete solution of the problem. 

Dividing the slab into zones, the rate equation R21 = (eq. 
I12> produces a similar expression for the i-th zone. 



5" + r]'p'S' = B{T' 



(26) 



with ry* and corresponding to the physical conditions in the zone. 
Inserting the expression for from eauation ll7l the CEP set of 
equations for the unknown 5" is 



(27) 



Since the factors /3' and M'-* depend only on optical depth, they are 
independent of the unknown variables (the zone source functions 
5") in this case. Therefore, the CEP technique transforms the two- 
level problem to a set of linear equations. This is a reflection of the 
linear relation between and S in eq.Qthat is maintained when 
the complete problem is handled in terms of optical depth as the 
independent variable. The CEP formulation produces directly the 
explicit linear equations in this case. 

It is interesting to note that the calculation of the mean intensity has also 
been done using Monte Cai'lo techniques (see, e.g., van Zadelhoff et al 2002 
and references therein). 
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Zones 



Time 



% error 



Semi-infinite 





SCP 


CEP 


SCP 


CEP 


20 


0.39 




36.3 


103.6 


40 


1.10 




23.9 


45.7 


100 


4.39 


.006 


10.9 


14.0 


200 


11.6 


.089 


5.5 


5.4 


600 


44.9 


1.70 


1.6 


1.2 



Table 1. Runtime (in seconds, on the same computer) required by the SCP 
and CEP methods to solve an atmosphere with e = 10~^ for the number 
of zones listed in the first column; in the SCP method this corresponds to 
the number of grid points. Omitted entries were too short for meaningful 
timing. The Hsted eiTor is the percent deviation from the result of an SCP 
calculation with 3,000 zones. 



We proceed now with solutions for semi-infinite atmospheres 
and finite-thickness slabs with constant physical conditions. When 
the temperature is constant, B{T) merely sets the intensity scale 
and only the dependence on e need be studied. 



3.1 Semi-infinite Atmosphere 

We start with the classical problem of a stellar atmosphere, where 
r measures distance from the surface and rt oo. The source 
function is subject in this case to the exact limits 



fe when r 







B X 



(28) 



1 when r 2> 1/e 



(e.g. Avrett & Hummer 1965). In order to capture both limit be- 
haviors we model the atmosphere as a slab divided logarithmically 
into 2 zones that cover ten orders of magnitude in optical depth 
from r — 10"'^ to rt — 10^, with the latter serving as a proxy 
for the atmospheric interior. The two faces of the slab are a-priori 
identical. When the radiative transfer equation is part of the cal- 
culation, this two-sided symmetry is broken by the boundary con- 
dition Iv ij ~ Tt,jj,) — 0, which introduces a radiation sink at 
the Tt-boundary. This is the case in ALI methods, including SCP. 
The CEP method, on the other hand, does not involve the radiation 
at all and thus cannot invoke boundary conditions to differentiate 
between the slab two faces. Instead, this is accomplished by the 
logarithmic division that starts at one end, and the great disparity 
that this introduces between photon escape from the two sides. The 
semi-infinite atmosphere could also be mimicked by doubling the 
slab with its mirror image and considering the source function only 
between one surface and the mid-plane. We have verified that the 
results of calculations with the two approaches are practically iden- 
tical. In order to compare the CEP method with SCP under identical 
conditions we present the results for logarithmic divisions increas- 
ing toward the slab surface at rt . 

Figure |3| shows the results for some representative models, 
ranging from e = 10"^ (TV = IQ-^A^cr) to e = 0.5 (iV = N^,). By 
example, the Call H line can be modeled in a 5000 K atmosphere 
with e = 3.65 xlO~^ (e.g., Avrett & Loeser 1987 and references 
therein). The top panel of each plot shows the solution obtained 
with the SCP method with 3,000 zones, displaying the proper limit 
behavior at both ends of the optical depth axis. The CEP method 
attains these solutions with a sufficient number of zones, validating 
our new technique. But the convergence with z is quite different for 
the two methods. 

At the surface, the SCP method is close to the exact solution 
already at z = 20 (only two zones per logarithmic decade) in all 
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Figure 3. The two-level model with various values of e (eg. 1241 in a semi- 
infinite atmosphere. The top panel of each plot shows the variation of the 
source function with depth into the atmosphere. The two other panels show 
the convergence to the exact solution as the number of zones z is increased 
for the CEP and SCP methods. Note the change in scale of the vertical axis 
between different panels. 
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cases; the deviation is less than 40% at e = 10"^ and decreases 
further as e increases. Increasing z brings rapid convergence. In 
contrast, deep inside, the rate of convergence is much more mod- 
erate. Furthermore, when e increases, both the magnitude of devi- 
ations and the rate of convergence around r ~ 10 remain almost 
the same for all e ^ 0.1. The behavior at both ends reflects the 
short characteristics nature of the method, in which only nearby re- 
gions are coupled, and the fact that the radiative transfer equation 
is always solved. Radiative transfer effects are minimal at small r, 
which is why the method attains easily the exact solution near the 
surface. But the effects are significant at the optical depths where 
the transition to thermalization occurs, the radiative transfer equa- 
tion must be repeatedly solved and the convergence in these regions 
is hardly improved by the increase in collision rates as long as 
remains sub critical. 

In an almost mirror behavior, at the surface the CEP method 
deviates from the exact solution by more than factor 2 at 2 < 100 
when e is small, and its convergence rate to the exact solution is 
slow there. However, deep inside the atmosphere, the deviations are 
actually smaller than at the surface. Moreover, when e increases, 
the deviations decrease everywhere. At e = 0. 1, the CEP method 
is within 7% of the exact solution everywhere already at 2 = 20; 
in contrast, this accuracy is attained by the SCP method only at 
z = 100. These properties are readily understood from the CEP 
formalism. Since the level population equations couple the entire 
atmosphere, the surface layers are affected by the behavior deep 
inside. And because the radiative transfer equation is avoided alto- 
gether, the method takes full advantage of the thermalization that 
approaches the surface when e increases. 

Performance statistics for the two methods are summarized in 
TableQfor the case e = 10^'^; the statistics for other cases show 
similar trends. The CEP technique attains the solution much faster 
than the SCP method in all cases. 

3.2 Finite-thickness Slabs 

One of the main coolants of photodissociation regions (PDR) is 
the 158 /im line of CII, whose emission is often modeled with a 
simple escape probability approximation of the two-level system 
(e.g., Tielens and HoUenbach 1985). When this approach employs 
the Capriotti expression for the escape probability (eq. I16> . it is 
identical to a CEP calculation with only one zone. For comparison 
with exact solutions, we include such single-zone CEP calculations 
in the results presented here. The numerical calculations employ z 
zones of equal thickness. 

Figure|4|shows the variation of the source function inside slabs 
of various optical thickness for e — 10~^. The displayed behavior 
is representative of all N <C N'^-^ cases. The results of single-zone 
CEP calculations provide reasonable approximations at small rt, 
but become poorer as the variation range of the source function 
gets wider with increasing rt. However, with only 20 zones the 
CEP results are within 1% of the exact solution everywhere when 
Tt ^ 10, 4% when Tt ^ 50 and 10% when Tt ^ 100. An accuracy 
better than 10% is always achieved when the optical thickness of 
each zone is < 5. In contrast, the SCP method does not reach this 
level of accuracy near the surface of a rt = 500 slab even with 200 
zones; as a grid- rather than zone-based method, it attempts to re- 
solve the surface structure even when that is not required. Further- 
more, SCP calculations require a large number of divisions even at 
moderate optical thickness; when rt = 10, a 10% accuracy requires 
100 zones. The reason, as noted above, is that the equation of ra- 
diative transfer must be solved repeatedly; the approach to thermal 
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Figure 4. The two-level model with e = 10^^ in slabs with overall optical 
thickness rt . The top panel of each plot shows the variation of the source 
function from the surface to the slab mid-plane in the exact solution and in 
a single-zone CEP calculation. The two other panels show the convergence 
to the exact solution with the number of zones z for the CEP and SCP 
methods. 



equilibrium of level populations deep inside the slab does not alle- 
viate this need, and large optical depths dictate a large number of 
zones. 

Since the CEP technique employs only level populations it 
takes full advantage of level thermalization. The difference with 
standard methods in the case of e = 0.5 (A'^ = N^^), shown in fig- 
ure |5| is striking. Already with one zone, CEP calculations pro- 
duce acceptable results inside every slab, even with rt as large as 
500; 20 zones suffice for 3% accuracy everywhere. In contrast, to 
achieve 10% accuracy, SCP requires 100 zones at a moderate rt = 
50, and even 200 zones are insufficient when rt = 500. The zone 
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Figure 5. Same as figure|4] only e = 0.5 



thicknesses in this case is r = 2.5, enough to challenge numerical 
solutions of the radiative transfer equation that SCP must perform. 

Table |2| summarizes the performance statistics in one repre- 
sentative case. The CEP method outperforms SCP by an even larger 
margin than in the case of an atmosphere. 



3.2.1 Slab emissivity 

The advantages offered by the CEP method are even more pro- 
nounced for slab emission calculations. Figure |S| shows a striking 
result: When e = 10"^ (i.e., 77 = 10^), the CEP method produces 
with a single zone the correct line cooling coefficient for all slabs 
with Tt ^ 500! This result is easy to understand. From the behav- 
ior of the function /3 at small and large r (Capriotti 1965) it follows 
that p(r) ~ 1 when r < 1 and that pij) ~ 1/r when r ^ 1, 
so that p(r) > 1 /rt . Therefore, slabs with rj > Tt have rjp > 1 
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Figure 6. The line cooling coefficient of a slab as a function of its overall 
optical depth in two-level models with various e. The top panel of each plot 
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The two other panels show the convergence to the exact solution with the 
number of zones z for the CEP and SCP methods. 
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Table 2. Performance comparison of the SCP and CEP methods, similar 
to Table in for a slab with e = 10~^ and Tt = 500. The listed errors in- 
clude the percent deviations of both the source function and the line cooling 
coefficient from the exact results. 



everywhere. From eQuations l25l andl5lit follows that in this case 
rt <r;: (29) 



VP 



-Bdr 

V 



The expression for j involves only input properties. That is, the line 
cooling coefficient can be calculated in this regime without even 
solving the problem. This result does not seem to have been recog- 
nized in the published literature. When the physical conditions are 
constant, j = Brt/ri; the source is optically thick yet its emission 
increases linearly with optical depth. The slab remains "effectively- 
thin" at large optical depths as long as rt ^ (i.e., ert <^ 1). And 
because the CEP method employs discretized forms of these ex- 
pressions, it reproduces the correct line emission irrespective of the 
division into zones. 

Since j{Tt) can be calculated without solving any equations, 
the single-zone calculation produces the correct emission even 
though it does not reproduce the correct population distribution — 
as is evident from both eq. 1291 and figure |4] the source function 
varies in the slab while the one-zone result is constant. Still, this 
constant value is just the right average to reproduce correctly the 
slab luminosity. Another perspective on this result is provided by 
the spectral shape of the emergent radiation, shown in figure^ The 
exact solution properly displays a self-absorption dip around line 
center (see, e.g., Avrett & Hummer 1965). The single-zone calcu- 
lation is incapable of producing this feature, but its flat-top shape 
does enclose the same area, reproducing the correct line luminosity. 
The simple one-zone calculation properly reproduces the overall 
number of photons emitted in the line, though not the frequencies 
where they emerge. 

When the problem is formulated in terms of r, eq. l29| gives the 
line emission directly from the input properties. When the problem 
is formulated instead in terms of densities and distances, eq. 1291 
implies that^ 



^ This result was noticed in the limit in which n2 <S ni by D. 
Neufeld in benchmark testing of radiative transfer codes, posted at 
I http://www.mpifr-bonn.mpg.de/staff/fvandertak/H20/radxf rtest.pdfi Note 
that the line cooling always obeys 



A = B21 J (C12 A^i - C2iN2)di 

= E21 j giCi2 (m - nae-^^i/''-'^) d£ 
as is evident from eqs. IslUlandlTol 
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Figure 7. Spectral shape of the flux emerging from slabs with e = 10~^ and 
optical depths as marked. The thick Unes are the result of the exact solution, 
the thin lines of a single zone CEP calculation. 



A = £'21 / giCi2 (ni 



712 



d£. 



(30) 



Although the condition r]p > 1 ensures that n2 <?C ni when 
i52i > fcT, n2 need not be negligible when £21 < kT. Therefore, 
the solution must be executed in this case to determine the popu- 
lation distribution and the actual value of Tt. Since the single-zone 
calculation does not produce the correct population distribution, its 
result for the overall optical depth can be wrong. To ensure the cor- 
rect assignment of Tt to the prescribed input, the problem must be 
solved properly, including the division into zones. 

When e increases, the slab ceases to be "effectively thin" and 
the line luminosity begins to deviate from the one-zone CEP result, 
as is evident from figure|5| Eventually, line thermalization sets in 
with further increase in e, and the single-zone result again becomes 
adequate. This behavior is further illustrated in figure|8| At a fixed 
Tt, the deviation from the single-zone CEP calculation reaches a 
maximum when e ~ 5/rt. When Tt = 50 the maximal deviation 
is ~ 20% at e ~ 0.1, when n = 500 it is ~ 70% at e ~ 0.01. 
Varying e away from that peak in either direction, the one-zone 
CEP calculation gives a progressively better approximation. 

With iV^, = 4x10^ cm"^ at T = 75 K, the 158 /xm line of 
ClI is in the regime > A^^r (e 0.5) in most cases of interest. 
Therefore single-zone CEP calculations for this line are expected 
to produce cooling rates accurate to better than ~ 10% under most 
circumstances. 



4 MULTI-LEVEL SYSTEMS 

Consider L energy levels. A trivial change from the two-level case 
is the addition of some bookkeeping indices. We designate level 
numbers with subscripts and zone numbers with superscripts. In 
zone i, the population per sub-state of level k is n], and the overall 
population is 



L 

9knl 

k = l 



(31) 



where gk is the level degeneracy. Unlike the two-level system, lo- 
cations in the slab cannot be specified by optical depth anymore 
because each transition has a different optical depth, which can be 
determined only after the unknown nj. are calculated. Instead, the 
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Figure 8. Deviation of tlie one-zone CEP result from the exact value of tlie 
line cooling coefficient as a function of e in two-level models with various 
slab thicknesses. 



partition into zones is done in terms of the geometrical distance 
from one surface. Denote by t the width of the i-th zone, then its 
optical thickness in the transition between lower level I and upper 
level u is 



(32) 



where E^i is the energy separation between levels u and I. The 
equivalent of eguation fTTl is then 



z,j \ k,k — l 



(33) 



when i > j. In complete analogy with eauations i 1 2111 7l and l 1 81 the 
level population equations are 



IT 



1=1 



AkiPkint + Cki luk - nie I 



(34) 



+ y 

^ gk 

u=fe + l 



AukPukn-u + {riu 



) 



Here 



Pui — Pui H i i_i 



Ml; 



(35) 
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Figure 9. Variation of the emission in the "^P cooling lines of OI with the 
number of zones in CEP calculations. Shown are the ratios to the exact 
solutions for slabs with T = 100 K, various H densities, as marked, and 
three representative oxygen column densities. At A^(0) 10^'' cm^^, the 
single-zone calculation produces the exact result for both lines. The critical 
density for each line is hsted in the top figure. 



-fa. 



(36) 



and where fi\i = P{t^i ^) and a'^f = t^^i I3(t'^i). This provides 
a set of L — 1 independent equations for the L unknown popula- 
tions in each zone, n\. Equation 13 ll for the overall density in the 
zone closes the system. The overall system of non-linear algebraic 
equations for the level populations in all zones is readily solved 
with the Newton method. 

It is convenient to switch to the scaled quantities n^/n* as the 
unknown variables and introduce the overall column density 



(37) 



Neither densities nor physical dimensions need then be specified 
since only J\f enters as an independent variable; the zone partition 
is done in terms of A/" rather than £. The problem is fully specified 
by three input parameters: density A'^ of collision partners and tem- 
perature T, which together determine the collision terms, and A/" 
(in fact, NIAvd), which sets the scale for all optical depths. 



4.1 Example — the O I cooling lines 

Together with the C II 158 jim, the lines of O I at 63 jim and 145 
/im dominate the gas cooling of warm PDR. Ratios and peak inten- 
sities of these lines are used to measure the gas density and temper- 
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Figure 10. The ratio r = j(145^m)/j(63/im) of the 01 coohng hnes as a function of oxygen column density for various temperatures and H densities, as 
marked. The top panel in each case shows the exact ratio r, the bottom panel the ratio of the results of single-zone calculations to the exact ones. 



ature (Tielens & Hollenbach 1985). In il3.2.1l we found that simple 
escape probability calculations do reproduce the proper Cn 158 

emission. We now examine the behavior of 01 lines through 
an exact CEP calculation of the three levels of the '^P system. We 
solve for slabs with constant physical conditions, specified by tem- 
perature, hydrogen density and oxygen column density M{0). 

Figure |9| shows the effect of varying the number of zones on 
the cooling lines emission at T = 100 K; the results are similar for 
T = 300 K and 500 K.*^ Single-zone CEP calculations produce the 
exact result at JV{0) ^ 10^^ cm~^, but deviate from it at larger 
column densities by amounts that increase with A/'(0). The devia- 

^ It is interesting to note that the 145 fj,m transition undergoes population 
inversion in the optically thin regime at low densities for temperatures above 
300 K. The reason is that the radiative lifetime of its upper level is more than 
five times longer than for its lower level. 



tion is different for each line, reflecting their different critical densi- 
ties, which are listed in the figure, and optical depths; for reference, 
at AfiO) = 10^^ cm"^ T(63^im) ~ 100 while r(145/im) varies 
from ~0.6 to ~2, depending on the density. Because of the differ- 
ent trends displayed by the emission in the two cooling lines, the 
one-zone calculation generally misses the exact result for their ratio 
r = j{145fim)/j{63fim). Fieure fTol shows the variation of r with 
column density and the deviations of the results of a one-zone cal- 
culation from the exact values. The single-zone results differ from 
the exact values by amounts that vary with temperature, density and 
oxygen column density. The deviations are generally largest at hy- 
drogen densities around 10^ - 10* cm~^. As is evident from this 
figure, the variation range of r is comparable to the error that can 
be introduced by its calculation with a single-zone. A reliance on 
such calculations can lead to erroneous conclusions regarding the 
physical conditions in a source. Indeed, from the observed ratio of 
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the two [01] lines and escape probability calculations, Caux et al 
(1999) deduced a mean gas temperature of 26±0.5 K, an H2 den- 
sity > 3x 10"' cm~^ and an [01] column density > 5x 10^^ cm~^. 
However, as is evident from fieure fTol the large values which they 
observed (r ~ 0.4) can be reached at lower columns for a wide 
range of somewhat higher temperatures. 

One-zone CEP calculations of the O I cooling lines are not as 
reliable as they are for the CII 158 /im. However, it does not take 
too many zones for the CEP method to achieve adequate accuracy. 
As is evident from figure|9| 20 zones suffice for accomplishing bet- 
ter than 10% accuracy, and the exact solution is reached to within 
1% with 40 zones in all cases. 



5 DISCUSSION 

The test cases presented here show that our new method not only 
provides an exact solution of the line transfer problem, it also out- 
performs the leading ALI solver by a margin even larger than that 
among different implementations of the ALI technique. Two fun- 
damental properties give the CEP method intrinsic advantages. The 
first is the calculation of J, the only radiative quantity required for 
solving the level populations. In the standard approach, I^i^j) is 
determined from a solution of eq.Qand J is calculated from eq.|3| 
in angular and frequency integrations that involve a-priori unknown 
dependencies on these two variables. Determining the dependence 
of Iv{ii) on V and /i is a major task for the solution of the radiative 
transfer equation. Deviations of the computed Iv{ii) from its exact 
angular shape and frequency profile contribute to the error in the 
computed J in each iteration. In contrast, the CEP method deter- 
mines J from the integration in eq.Qthat involves known depen- 
dencies on both u and /i; the dependence on these variables enters 
only from the optical depth r^(/i), and it is a-priori known from 
the input to the problem. The angular and frequency integrations 
are exact in the CEP method; the fact that the dependence of Iv{^i) 
on V and /i is unknown is altogether irrelevant. 

The other intrinsic strength of the CEP method is that it in- 
volves only level populations and thus takes full advantage of ther- 
malization wherever that sets in. In contrast, the ALI technique 
must repeatedly solve the radiative transfer equation in the entire 
source, even in thermalized regions, to determine the radiation field 
everywhere. 

5.1 Technicalities 

The great efficiency of the Newton method in solving non-linear 
equations is another advantage of the CEP method. The prerequi- 
site for a successful solution is a reasonable initial guess. An effi- 
cient strategy for working implementations of CEP is to start from 
the actual populations of a previous solution for similar physical 
conditions. A particularly useful approach is to start from the op- 
tically thin limit in which p = \ everywhere and solve the level 
populations from the corresponding linear equations. The column 
density N is now increased in small steps until the desired value is 
reached, using in each step the populations from the previous one 
as the initial distribution. This technique can also work in the oppo- 
site direction — start from thermal equilibrium populations and a 
very large N, and decrease N in small steps. An added advantage 
of this approach is that each step also provides information about 
the number of zones required for CEP convergence. 

The Newton method requires inversion of the Jacobian ma- 
trix, and the number of operations in this process increases as the 



third power of the matrix dimension. Although this rapid rise did 
not have a serious effect in the examples presented here, it could 
degrade the performance in cases of very large numbers of levels 
and zones. Matrix inversion is avoided in the iterative scheme de- 
signed by van der Vorst (1992) for solution of the linear system 
of the Newton method. In this scheme, geared toward sparse Ja- 
cobian matrices, only the non-zero matrix elements are stored and 
used. We have experimented with this method and found it to be 
quite useful for the CEP technique. It is particularly suitable for 
multi-level problems because they tend to produce sparse matrices, 
as each level generally couples to only a limited number of other 
levels. Other alternatives are to use quasi-Newton schemes (e.g., 
Broyden's method) like those employed by Koesterke et al. (1992), 
or evolve the set of differential equations 1341 until reaching steady 
state. 

The efficiency of CEP computations can be further enhanced 
with better grid design. Our solutions of the semi-infinite atmo- 
sphere employed grids with equal logarithmic spacing in r over ten 
orders of magnitude, resulting in extremely thick zones deep inside 
the atmosphere. For example, even with 2 = 600, the zone thickness 
was r = 1.7x10^ in the 10^ ^ r s; 10^ region. These extremely 
thick zones do not pose any difficulties to CEP computations be- 
cause they occur in regions where the populations are thermalized. 
Indeed, the zones could be even thicker in a much larger fraction of 
the source without compromising accuracy. It should thus be pos- 
sible to achieve the same accuracy with fewer zones by concentrat- 
ing them in the regions were the populations deviate from thermal 
equilibrium. Since the number of zones is the single most impor- 
tant factor in determining CEP runtime, a more sophisticated grid 
construction will make the method even more efficient. We intend 
to investigate the implementation of adaptive gridding algorithms 
in future work. 

An additional increase in efficiency can be easily gained in 
practical applications that do not require the extreme precision we 
imposed in this comparative study. Here the functions a and /3 were 
calculated using the integral definition in eq. 1221 repeatedly per- 
forming highly accurate quadrature. Instead, one could employ the 
approximate series expansion derived by Capriotti (1965) for the 
function 13. We have verified that this rapidly convergent series is 
always within 3% of the exact result.^ Another option is to calculate 
once a finely spaced table of the a and /3 functions and interpolate 
between its elements with an efficient algorithm. 

5.2 CEP and ALI 

The CEP method is suitable for solution also with the ALI ap- 
proach. Starting from eq.|3|in the operator form J = AS, the ALI 
technique is based on the operator splitting A = A* + (A — A*), 
where A* is an approximation to the A operator. The mean inten- 
sity is obtained from the approximate expression J = A*5'+(A — 
^(.■j^prcv jjj^j involves the source function from the current and 
previous iterations. This approximation becomes exact upon con- 
vergence, when S — S'^'^°^ . As already noted, it has been shown 
that the optimal choice for A* is the diagonal of A. 

Equation |6| gives J = (1 — p)S, the A-operator form of the 
CEP method. From eq. 1171 the matrix elements of the CEP A- 
operator are simply 

^ In Capriotti (1965), the small- and large-r portions of the expansion were 
joined at line-center optical depth tq = 5. They should be joined instead at 
TO = 3.41 for a smooth transition. 
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A,, = (l-/3')5.,-— ^(1-5,,) (38) 

where 5ij is the Kronecker delta. Thanks to the known dependence 
on V and fi in the CEP approach, this expression allows the usage 
of approximate operators A* of increasing complexity without any 
additional computational effort in the calculation of the matrix el- 
ements. In standard ALI techniques A* is the diagonal of A, i.e., 
A*j- = (1 — We have implemented this choice in an ALI 

solution of eq.|23|and performed the two-level model calculations 
presented in this paper also with this technique. In all cases, the so- 
lution converged to the exact same results as the algebraic CEP 
equations 1271 Runtime for this ALI implementation of the CEP 
method was on par with the SCP method up to 200 zones, but fell 
behind at larger z.^ Given that in the CEP ALI implementation we 
have adopted the standard choice for A*, the optimal choice in the 
CEP approach could well be different, improving the performance. 

It is also important to point out that the CEP method could, 
in principle, be implemented in the framework of the Gauss-Seidel 
and Successive Overrelaxation iterative methods. These methods 
were first applied in radiative transfer problems by Trujillo Bueno 
& Fabiani Bendicho (1995) using SCP as the formal solver. They 
can lead to an order of magnitude improvement in the number of 
iterations used to reach convergence, with a time per iteration that is 
virtually the same as the method based on the lacobi iteration. Also 
of interest is the possibility of implementing the CEP method in the 
linear (Steiner 1991) or the non-linear (Fabiani Bendicho, Trujillo 
Bueno & Auer 1997) multi-grid methods. All of these issues will 
be addressed in future investigations. 



5.3 Extensions 

All the examples presented in this paper involved constant physi- 
cal conditions. Variable conditions are handled by simply starting 
with zones that have constant physical conditions and proceeding 
to refine those divisions as required by the CEP solution accuracy. 
EQuation l34l for the level populations already incorporates the han- 
dling of variable conditions by allowing the temperature and colli- 
sion rates to vary between the zones. 

For simplicity, our method was introduced in the context of a 
quiescent slab with the Doppler line profile. None of these simpli- 
fications represents an inherent limitation of the CEP method. The 
formal expressions do not specify the shape of ^{x), and other line 
profiles can be implemented just the same. Extension from the slab 
to other geometries, although straightforward, requires some more 
work. Thanks to the planar symmetry, the angular variation of op- 
tical depth in a slab is simply r (/i) = t(/x = independent of 
either position or density profile. This symmetry does not carry to 
any other geometry; even in the case of spherical symmetry, the an- 
gular variation of t cannot be calculated at any point other than the 
center without specifying the density profile. However, generaliz- 
ing the basic CEP relation eq.0to handle any geometry is straight- 
forward, and the fundamental advantage of integration over known 
frequency and angular variations remains intact. Finally, recalling 
that large velocity gradients was the context in which the escape 
probability approximation was originally introduced by Sobolev, 
the CEP method is well suited for exact handling of this case too. 

The escape probability approach has been used in a number of 



° We implemented the Ng (1974) acceleration technique to improve the 
convergence rate of the ALI method in both SCP and CEP. 



simplified calculations of complex problems. These include: over- 
lapping of spectral lines (so called "line fluorescence"), important 
for various ionic transitions (e.g., Bowen lines) in photoionized 
regions and OH lines in molecular clouds (Guilloteau, Lucas & 
Omont 1981, Elitzur & Netzer 1985, Lockett & Elitzur 1989); the 
effect of line overlap with underlying continuum (Netzer, Elitzur & 
Ferland 1985); and photoionization (Elitzur 1984). Importing these 
applications into the CEP framework is straightforward. Finally, 
another extension is the application of the CEP method to the self- 
consistent solution of the radiative transfer equations for polarized 
radiation and of the statistical equilibrium equations for the density 
matrix, the so-called non-LTE problem of the 2"'^ kind (see, e.g., 
Landi DeglTnnocenti 2003). We plan to provide these extensions 
in future publications. 

5.4 Conclusions 

While our new method outperforms the current leading techniques, 
its greatest advantage is its simplicity and ease of implementation. 
The CEP method employs a set of algebraic equations (eq. l34> that 
are already incorporated in numerous widely used codes based on 
the escape probability approximation. All that is required for an 
exact solution with these existing codes is to augment the escape 
probability with the zone-coupling sum on the right-hand- side of 
eg. 1351 With this simple modification, the multi-level line transfer 
problem is solved exactly. 
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APPENDIX A: EXTERNAL RADIATION 

The only effect of external radiation on the rate equations is to mod- 
ify the exchange rate ii^j between levels u and I in the i-th zone 
(see eQ. ll2> according to 

K,^ Ki~ B„iJl(n'^~nl), (Al) 

where Jl is the zone average (as in eQ. ll3> of the contribution of 
the external radiation to the local J. When the external radiation 
corresponds to the emission from dust which permeates the source, 
Jl is simply the angle-averaged intensity of the local dust emission 
in the i-th zone. When the external radiation originates from outside 
the slab and has an isotropic distribution with intensity Is (= Je) 
in contact with the r = face, then 

J^ = hJ^^i<?-<i''°)- (A2) 

When the slab is illuminated by parallel rays with intensity Je (= 
4:11 J e) entering at direction (/io, 0o) to the r = face then 

Jl = -^e [TKi/Mo) - 7(Tu7Vmo)] (A3) 

'''ul 

where 

[l _ dx. (A4) 

- OO 
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